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ABSTRACT 

The magnetorotational instability (MRI) drives vigorous turbulence in a region of protoplanetary disks where 
the ionization fraction is sufficiently high. It has recently been shown that the electric field induced by the 
MRI can heat up electrons and thereby affect the ionization balance in the gas. In particular, in a disk where 
abundant dust grains are present, the electron heating causes a reduction of the electron abundance, thereby 
preventing further growth of the MRI. By using the nonlinear Ohm’s law that takes into account electron 
heating, we investigate where in protoplanetary disks this negative feedback between the MRI and ionization 
chemistry becomes important. We find that the “e-heating zone,” the region where the electron heating limits 
the saturation of the MRI, extends out up to 80 AU in the minimum-mass solar nebula with abundant submicron¬ 
sized grains. This region is considerably larger than the conventional dead zone whose radial extent is ~ 20 AU 
in the same disk model. Scaling arguments show that the MRI turbulence in the e-heating zone should have 
a significantly lower saturation level. Submicron-sized grains in the e-heating zone are so negatively charged 
that their collisional growth is unlikely to occur. Our present model neglects ambipolar and Hall diffusion, but 
our estimate shows that ambipolar diffusion would also affect the MRI in the e-heating zone. 

Keywords: accretion, accretion disks - instabilities - magnetohydrodynamics (MHD) - planets and satellites: 
formation - plasmas - protoplanetary disks - turbulence 


1. INTRODUCTION 

Magnetorotational instability (MRI; Balbus & Hawley 
1991) is widely regarded as a mechanism driving turbulence 
in protoplanetary disks. Vigorous MRI turbulence provides an 
effective viscosity that allows disk accretion at a rate consis¬ 
tent with observations (Hawley et al. 1995; Elock et al. 2011). 
MRI also generates outflows from the disk surface (Suzuki & 
Inutsuka2009; Suzuki et al. 2010; Bai 2013; Lesur et al. 2013; 
Eromang et al. 2013). In addition, MRI turbulence have many 
important effects on the evolution of solid particles and planet 
formation. The effects include diffusion of small grains (Car- 
ballido et al. 2005), concentration of larger solid particles (Jo¬ 
hansen et al. 2006), enhancement of particles’ relative veloc¬ 
ity that could lead to their collisional disruption (Carballido 
et al. 2010) of meter-sized bodies, and random migration of 
planetesimals (e.g., Laughlin et al. 2004; Nelson & Gressel 
2010). 

However, in weakly ionized protoplanetary disks, the satu¬ 
ration level of MRI turbulence depends strongly on non-ideal 
MHD effects and hence on the ionization state of the disks. 
Since thermal ionization is relevant only close to the cen¬ 
tral star (Umebayashi 1983), the dominant part of the disks 
is ionized only by high-energy sources such as stellar X-rays 
(Glassgold et al. 1997) and galactic cosmic rays (Umebayashi 
& Nakano 1981). Deep inside the disks, the ionization frac¬ 
tion is significantly low because these ionizing radiations are 
attenuated and because recombination proceeds fast. 

The low ionization fraction gives rise to fast Ohmic dissipa¬ 
tion that stabilizes the MRI (Sano & Miyama 1999). Such a 
region is called the “dead zone” (Gammie 1996). The MRI 
is also suppressed by ambipolar diffusion near the surface 
of the disks (Desch 2004; Bai & Stone 2011; Dzyurkevich 
et al. 2013). The Hall effect can either stabilize or destabi¬ 


lize the MRI depending on the orientation of the magnetic 
field relative to the disk rotation axis (Wardle 1999; Wardle & 
Salmeron 2012; Bai 2014). 

A number of studies have quantified how far the dead zone 
extends in protoplanetary disks. Gammie (1996) assumed that 
the MRI is stable in a region where the column density ex¬ 
ceeds the attenuation depth 100 g cm“^) of galactic cos¬ 
mic rays. More sophisticated models that incorporate ion¬ 
ization and recombination (e.g., Sano et al. 2000; Semenov 
et al. 2004; Ilgner & Nelson 2006; Bai & Goodman 2009; 
Dzyurkevich et al. 2013) showed that the MRI can be inac¬ 
tive even at lower column densities, with the predicted dead 
zone extending to ~ 20 AU from the central star when micron¬ 
sized dust grains are abundant in the disks. The abundance of 
small grains is relevant here because these particles efficiently 
sweep up plasma particles and thus lower the ionization frac¬ 
tion. 

All the previous studies on the dead zone assumed that vig¬ 
orous MRI turbulence is sustained outside the zone. However, 
Okuzumi & Inutsuka (2015, henceforth 0115) suggested that 
the ionization fraction would be decreased by electric fields 
induced by MRI turbulence. In a magnetorotationaly unstable 
region, the MRI turbulence generates strong electric fields as¬ 
sociated with the growth of magnetic fields. Plasma particles 
are accelerated by the strong electric fields and are scattered 
isotropically by collisions with neutral gas particles, leading 
to increase of their thermal velocity. In particular, electrons 
are more easily heated compared to ions because light parti¬ 
cles are easily scattered. Therefore, the sufficiently developed 
electric fields of MRI turbulence increase electron tempera¬ 
ture in a weakly ionized gas {electron heatings Inutsuka & 
Sano 2005). The heated electrons frequently collide with and 
stick to dust grains. As a result, the electron heating decreases 
the ionization fraction. 
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Reduction of ionization fraction caused by the electron 
heating amplifies Ohmic dissipation, and, as a result, MRI tur¬ 
bulence may be suppressed. This negative feedback that the 
MRI growth causes suppression of the MRI can be a saturat¬ 
ing mechanism of MRI. Although previous simulations (e.g., 
Sano & Stone 2002; Simon et al. 2015; Flock et al. 2015; 
Bai 2015) in a well ionized regions showed that MRI turbu¬ 
lence sustains a fully developed state, the turbulence strength 
may be suppressed at a lower saturation level by the effect of 
the electron heating. However, it is not clear whether elec¬ 
tric fields can sufficiently grow to heat up electrons before the 
MRI fully develops, and whether the decrease of the turbu¬ 
lent saturation level is meaningful. The goal in this paper is 
to investigate where in protoplanetary disks the electron heat¬ 
ing affects MRI turbulence and estimate how the saturation 
level would be suppressed. This investigation is the first-step 
towards exploring the importance of the electron heating in 
protoplanetary disks. In this study, we take into account only 
the Ohmic dissipation and neglect the other non-ideal effect 
of MHD for simplicity. 

In Section 2, we present the disk model, simplified plasma 
heating model, and ionization balance. In Section 3, we 
present some conditions for MRI growth and some criteria for 
mapping of turbulent state in a disk. We also briefly summa¬ 
rize the turbulent state and calculation steps. In Section 4, we 
show where the electron heating affects MRI turbulence. We 
also consider cases with various parameters. In Section 5, we 
estimate how the electron heating suppresses MRI turbulence. 
In Section 6, we discuss the effect of heated electrons on the 
electric repulsion and the collisional growth of dust grains. In 
Section 7, we discuss the effects neglected in our study. In 
Section 8, we present a summary of the results. 

2. DISK AND IONIZATION MODELS 
2.1. Disk Model 

We consider a gas disk around a solar-mass star. We assume 
that the surface density of the disk gas obeys a power law 

E(r) = 1.7 X IOVe (y^) gcm-2, (1) 

where r is the distance from the central star, and fz is a di¬ 
mensionless parameter. The choice offz = l corresponds to 
the minimum-mass solar nebula (MMSN) model of Hayashi 
(1981), which we take as the fiducial model. 

We assume that the disk is optically thin and give the tem¬ 
perature profile as (Hayashi 1981) 


/ r \-i/2 

where the central star is assumed to have the solar luminosity. 

The sound speed is given by Cg = ^^kTJm^, where is the 
mass of a neutral gas particle, and k is the Boltzmann constant. 
Assuming = 2.34 amu and using Equation (2), we have 


where pc is the mid-plane density and H = is the gas 
scale height with Q = 2.0 x 10“^ (r/1 AU)“^/^ s“^ being the 
orbital frequency (note that a solar-mass star is assumed). Us¬ 
ing the relation E = f_^pdz = ^^J^Hpc, we have 

Pcif) = 1.4 X 10"Ve (y^) 8 

Thus, the number density of gas particles tin = plnin is given 
as 

n„(r,z) = 3.5x IO'VsX^Y;^) 


As we will describe in Section 3.1, the criteria for MRI de¬ 
pends on the magnetic field strength in the disk. Following 
Sano et al. (2000), we consider a net (large-scale) vertical 
field 5^0 threading the disk and specify its strength with the 
plasma beta at the midplane, Pc = ^npcC^J If we use 
Equations (3) and (5), the net vertical field strength can be 
expressed as 


B,o(r) = 0.59fl^^ 



G. 


(7) 


For simplicity, we will assume that Pc is constant in the radial 
direction. 

The charge reaction model adopted in this study takes into 
account the effects of grain charging on the ionization bal¬ 
ance. For simplicity, we assume that dust grains are well 
mixed in the gas so that the dust-to-gas mass ratio fdg is a 
global constant. We also assume that the grains are spheri¬ 
cal and single-sized with radius a (taken as a free parameter) 
and internal density p, (fixed to be 3 g cm“^). From these as¬ 
sumptions, the number density of dust grains rid is given by 
3fdgpl{4na^p.), which is expressed as 


nd(r,z)=lAxWfz 
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The disk is assumed to be ionized by galactic cosmic rays, 
stellar X-rays, and radionuclides. The ionization rate can be 
expressed as 

f = fcR + fxR + ^RN, (9) 


where fcR. fxR, and ^rn stand for the contributions from 
cosmic rays. X-rays, and radioactive decay, respectively. 
The cosmic ray distribution is expressed as (Umebayashi & 
Nakano 2009) 


. fcR,o j / T 
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3/ r -t 

Cs(r) = 1.0 X 10^^Y^j cms ‘. (3) 

We assume that the gas disk is hydrostatic in the vertical 
direction and give the vertical distribution of the gas density 
as 

p(r, z) = Pc{r) exp j- (4) 


where fcR,o = 1.0 x 10 s Ms the characteristic ionization 
rate of cosmic rays, ;^(r, z) = f p(r, z')dz' is the vertical gas 

'Jz 

column density above height z, and xcr = 96 g cm“^ is the 
attenuation depth of ionizing cosmic rays. The ionization rate 
of X-rays is expressed as (Bai & Goodman 2009) 

^XR=i029tgs-i (tXu) 
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where ;^XR,i and;^XR ,2 are taken to be 6 x 10“^ g cm"^ and 
3 g cm“^ respectively, fxR,i and ^xr ,2 are taken to be 6 x 
10“^^ s“^ and 1 x 10“^^ s“^ respectively. We take = 
2 X 10^^ erg s“^ in accordance with the median X-ray lumi¬ 
nosity of solar-mass young stars (Wolk et al. 2005). The ion¬ 
ization rate of the radionuclide is expressed as (Umebayashi 
& Nakano 2009) 


(~ 10 eV). The results of our calculations show that this as¬ 
sumption holds in most parts of protoplanetary disks. 

We denote the mean drift velocity and mean kinetic en¬ 
ergy of a charged species a (= i for ions, e for electrons) 
by {Va) and respectively. In a weakly ionized gas with 
an applied electric field E, the momentum and energy of the 
charged species are determined by the balance between the 
neutral gas drag and acceleration by the electric field (Her- 
shey 1939). Explicitly, the solution of the momentum and 
energy balance equations can be written as (Equations (A9) 
and (AlO) of 0115) 

nify rfin 

(Va) = — - -qaEAta, (13) 

maMn 


^RN = 7.6xlO-'5|^j s-‘. (12) 



{mg + mnf 

2{mamnY 


iqaEAtg) , 


(14) 


2.2. Simplified Plasma Heating Model 

As we will describe in Section 3.1, the criterion for MRI 
depends on the ionization fraction in the disk. We employ 
a simple ionization model proposed by 0115 to calculate the 
ionization fraction taking into account plasma heating by a 
strong electric field. The model determines the ionization 
fraction of the gas at each location of a disk from the balance 
between ionization by external high-energy sources (e.g., cos¬ 
mic rays and X-rays), recombination in the gas phase, and ad¬ 
sorption of ionized gas particles onto dust grains. The rates of 
recombination and adsorption generally depend on the tem¬ 
peratures of ions and electrons, Ti and Tg. Previous ionization 
models assumed that Ti and Te are equal to the neutral gas 
temperature T. By contrast, the model of 0115 determines 
Ti and Tg as a function of the electric field strength E. Eor 
simplicity, positive ions are represented by the single species 
HCO^, which is good as a first-order approximation when 
heavy molecular ions that recombine through dissociation re¬ 
actions dominate (Umebayashi & Nakano 1990; Dzyurkevich 
et al. 2013). We do not consider negative ions. Although pro¬ 
duction of negative ions is rare in cool protoplanetary disks, 
electrons heated to > 3 eV can produce negative hydrogen 
ions H“ via dissociative electron attachment H 2 +e“ ^ H“ -rH 
(Wadehra 1984). However, H“ would be instantly destroyed 
by CO, the most abundant molecule after H 2 , via the reaction 
H“ -h CO ^ HCO -r e“ (Eerguson 1973). Eor this reason, we 
may safely neglect the dissociative electron attachment during 
electron heating. 

In this study, we make two further simplifications to the 
original model of 0115. Eirstly, we calculate the electron tem¬ 
perature Te by solving the equations of momentum and en¬ 
ergy conservation rather than by using the solution to the full 
Boltzmann equation. The rate coefficients for gas-phase re¬ 
combination and plasma adsorption onto grains are then eval¬ 
uated by approximating the velocity distribution function with 
a Maxwellian with temperature Tg. The approach greatly sim¬ 
plifies the analytic expressions of the rate coefficients that oth¬ 
erwise involve confluent hypergeometric functions (see Sec¬ 
tion 3 of 0115). Such an approach was originally proposed 
by Hershey (1939) for calculating the mobility of heavy ions 
at a high electric field, and 0115 followed this approach to 
compute the ion temperature T/. In this study, we apply this 
approach to both U and Tg. Secondly, we neglect the im¬ 
pact ionization of neutral molecules by electrically heated 
electrons by assuming that the electron energy in MRI tur¬ 
bulence is well below the ionization potential of the neutrals 


where and Atg are the charge, mass and mean free time 

of the plasma particles (e.g., qe = -e and qi = e, where e is 
the elementary charge). Since the magnetic field is neglected 
in this study, the mean drift velocity is parallel to the electric 
field. In a weakly ionized gas, the plasma mean free time is 
determined by neutrals gas particles, 

Atg = (jlfi {O'gnVgn}) , (15) 

where Vg^ is the relative velocity between a plasma particle 
and a neutral particle, and o-g^ is the momentum-transfer cross 
section for the plasma-neutral collision. Eor electrons, CTen is 
approximately constant at low energies (Yoon et al. 2008), 
and therefore we may approximate {o-enVen) as aen (ven)- For 
ions, {o-inVin) is approximately constant owing to the polariza¬ 
tion force between ions and neutrals (Wannier 1953). Equa¬ 
tions (13) and (14) are exact only when Atg is constant, but 
still hold in a good accuracy even when Atg is velocity- 
dependent (Wannier 1953). 

The plasma temperature Tg is defined so that ?>kTgl2 
is equal to the kinetic energy of random motion, {eg) - 
(t^a)^ /2. Using Equations (13) and (14), Tg can be written 
as 

rr _rr. + ^ 

Ea — T -\- (qgEAtg) . (lb) 

3km^mn 

/ 

Eor electrons, we approximate (Ven) in Ate with = 

sl^kTelme. This allows us to solve Equation (16) with respect 
to Te, and we obtain 


where 


Te = T 
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(18) 


is the critical field strength above which electron heating be¬ 
comes significant. We have assumed in deriving 

Equation (17). Eor ions. Equation (16) directly gives 


Ti = T 


= T 


1 -r 


2(mi + mn)^me cr^n^T 
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where we have set {(TinVin) = 1.6 x 10“^ cm^ s“^ (Nakano & 
Umebayashi 1986) and (Ten = 10“^^ cm^ (Yoon et al. 2008) in 
the second expression, and used m/ = 29 amu. 


2.3. Ionization Balance and Accuracy of Simplified 
Approach 

We calculate the plasma densities in a protoplanetary disk 
taking into account grain charging. The equations that de¬ 
scribe the ionization balance in a dusty disk are (Equations 
(32), (33) and (35) of 0115) 


^nn - K,^c(Te)nine - Kdeif, Te)ndne = 0 , ( 20 ) 

^nn - Kr^c{Te)nine - Kdi((p, Ti)ndni = 0, (21) 

ni - + Znd = 0, (22) 


where ne and ni are, respectively, the number density of elec¬ 
trons and positive ions; ^rec is the gas-phase recombination 
rate; Kde and Kdi are the adsorption rates of electrons and 
ions onto grains; Z is the grain charge number; and f is the 
coulomb potential on grain surface. 0 is related to Z as 

eZ 

0=—. (23) 

a 


As the collisional frequency, ^rec and Kde depend on the elec¬ 
tron temperature while Kdi depends on the ion tempera¬ 
ture Ti. Kde and Kdi also depend on the coulomb potential of 
a grain surface 0. For HCO^, the recombination rate ^rec is 
given by (Ganguli et al. 1988) 


Kr^cTe) = 2.4 X 10-’ 



3 -1 

cm s 


(24) 



Figure 1. Test of the simplified plasma heating model presented in Section 
2.3. The solid curve shows the J-E relation for ‘model C’ of 0115 derived 
using the exact electron velocity distribution (see Figure 10 of 0115), while 
the dashed curve shows our reproduction based on the simplified approach. 

the procedure presented by Okuzumi (2009, their Section 2.2; 
see also Section 3.2.4 of 0115). 

To test the accuracy of our simplified approach, we repro¬ 
duce the current-field relation including plasma heating (the 
nonlinear Ohm’s law of 0115) with adopting the calculation 
steps in 0115. Current density is generally given by 


J{E) = qene (Ve) + qini (Vi) . (27) 


Approximating the ion velocity distribution by a Maxwellian 
with mean velocity (u/) and temperature Ti, Kdi is given by 
(Shukla & Mamun 2002, 0115) 


Kdii(f>,Ti)=na^ 



2kTi ) 



kTi -f 2^101 
ttli {Vif 


erf 


^jTFTiJmi 


. (25) 


In this study, we also approximate the electron velocity dis¬ 
tribution by a Maxwellian with temperature Te. The drift ve¬ 
locity {Ve) can be neglected here since the drift speed \{Ve)\ 
is generally much smaller than the random speed ~ ^^F^Jn^e 
owing to the smallness of Melmn (see Golant et al. 1980; Lif- 
shitz & Pitaevskii 1981). The electron adsorption rate coeffi¬ 
cient Kde is given by the simple expression (Shukla & Mamun 
2002) 


Kdeitp, Te) = n(Z 



1 + ^1. ^>0, 



(26) 


Including plasma heating, the number densities depend on the 
electric fields strength E. To obtain the current density, we 
first calculate plasma temperatures Te and Ti from Equations 
(17) and (19) in an applied electric field E. We then calcu¬ 
late the number densities of plasma ne and ni from the ion¬ 
ization balance (Equation (22)). We finally obtain the cur¬ 
rent density using Equations (13) and (27). In Figure 1, we 
compare our result with the result of 0115 for the parameter 
set ‘model C’ of 0115. We find that our calculation reason¬ 
ably reproduces the previous result even at high field strengths 
{E > 10“^ esu cm“^) where electron heating is significant. 
The maximum relative difference between the two results is 
37%. 


3. ACTIVE, DEAD, AND E-HEATING ZONE 
3.1. Conditions for MRl Growth 

In the limit of ideal MHD, the criterion for the MRI is given 
by (Balbus & Hawley 1991) 

dideal < fi, (28) 

where 

Aideai = 2;r^ (29) 


It should be noted that Equations (26) and (25) assume per¬ 
fect sticking of ions and electrons onto grain surfaces. This 
is a good approximation as long as the plasma temperatures 
are well below 100 eV (see Section 3.2.2 of 0115 for more 
discussion). 

Equations (20)-(22) determine ne, ni and Z at each location 
in a disk as a function of E. We solve these equations using 


is the characteristic wavelength of the most unstable axisym- 
metric MRI modes, and vaz = and are the vertical 

components of the Alfven velocity and magnetic field, respec¬ 
tively. Equation (28) expresses that the MRI operates when 
the lengthscale of the MRI modes is smaller than the vertical 
extent of the disk. When viewed as a function of z, dideai in¬ 
creases with z because p decreases toward the disk surface. If 
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we use Equation (4), the above MRI criterion can be rewritten 
in terms of height as 

z < ^21n08J8;r2) H = Hweai, (30) 

where the height i^^ideai defines the upper boundary of the MRI 
active zone. 

Inclusion of a finite Ohmic resistivity rj introduces another 
criterion for MRI growth. The criterion can be expressed in 
terms of the Elsasser number (Turner et al. 2007) 


The instability grows when 

A > 1 (32) 

and decays when A < 1 (e.g., Sano & Miyama 1999). 

3.2. Zoning Criteria 

Here we describe how to determine turbulent state at a po¬ 
sition in protoplanetary disks. Electron heating affects on the 
MRI turbulence when the ionization fraction is sufficiently de¬ 
creased. We express the condition that the heating takes place 
and affect MRI turbulence, and then summarize three turbu¬ 
lent states of MRI and steps of zoning a disk into the state. 

Eor electron heating to take place, the field must be suffi¬ 
ciently amplified before MRI turbulence reaches a fully de¬ 
veloped state that means the stop of MRI growth. Muranushi 
et al. (2012) performed a local unstratified resistive MHD 
simulation and found that the fully developed current density 
is 



where /sat ~ 10 according to the results by Muranushi et al. 
(2012). Here, we assume /sat to be /sat = 10 and the maxi¬ 
mum current density is /max- Thus, when the current density 
reaches /max before electric field reaches the criterion for elec¬ 
tron heating MRI turbulence does not cause the electron 
heating. 

As we will describe later in this section, we use current 
density to decide whether electron heating take place or not. 
Therefore, we transform the condition for suppressing MRI 
into a form using current density. We adopt A = 1 (Equation 
(32)) as the criterion for suppressing MRI which is triggered 
by electron heating. Using the electric conductivity (Tc and 
the relation rf = jAnac, the condition for sustaining MRI 
turbulence A > 1 leads to a condition ac > c^Q.I{Anv\^. Un¬ 
der the Ohm’s law J{E) = acE, the condition can be rewritten 
as a lower limit to the current density 

J(E) > Ja=i(E\ (34) 

where ^ 

Ja=i(E) = (tAA = l)E = -^E. (35) 

Using the above criteria, we can classify a region in pro¬ 
toplanetary disks into three different zones corresponding to 
three turbulent state of MRI. 

1. Dead zone. Because of the low ionization fraction. 
Ohmic dissipation suppress all the unstable MRI mode. 
Suppressed MRI does not generate turbulence and also 



Figure 2. Flow chart showing key steps of zoning a protoplanetary disk into 
the dead, active, and e-heating zones. 


current density. We will refer to the region where 
MRI is completely suppressed as the “dead zone”. In 
this case, the condition of Ohmic dissipation (Equation 
(35)) is satisfied with no MRI turbulence. 

2. E-heating zone. Electric fields of MRI turbulence be¬ 
come sufficiently high for electron heating to be caused. 
The Ohmic dissipation is amplified by the electron heat¬ 
ing after the MRI grows. We will refer to the region 
where electron heating affects MRI turbulence as the 
“e-heating zone”, where the “e” refers to both “electric 
field” and “electron.” In this case, current density falls 
down the critical current density of Ohmic dissipation 
(Equation (35)). 

3. Active zone. MRI sustains fully developed turbulent 
state because the gas is sufficiently ionized so that 
Ohmic dissipation is not efficient. We will refer to the 
region where vigorous MRI turbulence is sustained as 
the “active zone” in this study. In this case, the current 
density / reaches and sustains its maximum value /max 
before electron heating reduces the MRI turbulence. 

We summarize the calculation steps for zoning the disk re¬ 
gion under some assumptions. We assume that the electric 
field strength correspond to the activity of MRI turbulence 
since developed MRI generates strong electric fields. The 
growth of MRI implies increasing electric fields, and the de¬ 
cay of MRI implies decreasing electric fields. Eurthermore, 
we also assume that magnetic fields are not varied by the MRI 
growth for simplicity. Under these assumptions, we deter¬ 
mine the turbulent state at the position with following steps 
(see Eigure 2): Eirst, we select a calculated position in the 
region satisfying Equations (30) of a disk. We then calcu¬ 
late values at the position with setting E = Q. When MRI 
is initially suppressed by Ohmic dissipation, i.e., A < 1 at 
£■ = 0, the positions belong to the dead zone. During satisfy¬ 
ing unstable condition, i.e., A > 1, the electric field strength 
E is increased from £" = 0 with iterating until the turbulent 
state at the position is determined. We calculate current den¬ 
sity J{E) and assess some conditions in E. When MRI tur¬ 
bulence causes electron heating and Ohmic dissipation be¬ 
come efficient, i.e., J{E) = Ja=i, the position belongs to the 
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Table 1 

Sizes of the Dead and E-heating Zones for Various Parameter Sets 


A 

a (jdm) 

fdg 

h 

Outer radius (AU) 

Dead zone E-heating zone 

10^ 

0.1 

10-2 

1 

18 

74 

10^ 

0.1 

10-2 

1 

24 

82 

10“* 

0.1 

10-2 

1 

34 

82 

10^ 

0.1 

10-2 

1 

56 

82 

10^ 

0.1 

10-2 

1 

24 

82 

10^ 

1 

10-2 

1 

11 

39 

10^ 

10 

10-2 

1 

8 

19 

10^ 

100 

10-2 

1 

8 

11 

10^ 

0.1 

10-1 

1 

52 

151 

10^ 

0.1 

10-2 

1 

24 

82 

10^ 

0.1 

10-3 

1 

12 

41 

10^ 

0.1 

10-4 

1 

8 

20 

10^ 

0.1 

10-2 

10 

55 

149 

10^ 

0.1 

10-2 

3 

36 

114 

10^ 

0.1 

10-2 

1 

24 

82 

10^ 

0.1 

10-2 

0.3 

14 

44 



Figure 3. Cross-section view of the fiducial protoplanetary disk indicating 
the location of the dead, e-heating and active zones (red, green-shaded, and 
blue regions, respectively). The dashed line shows the gas scale height //, 
while the dotted line shows the critical height //ideal below which the MRI 
criterion in the ideal MHD limit is satisfied (see Equation (30)). 

e-heating zone. When MRI turbulence is fully developed, i.e., 
J{E) = /max, the position belongs to the active zone. We con¬ 
duct the above steps in the whole region in a disk, and zone a 
protoplanetary disk into the dead, active, and e-heating zones. 

4. LOCATION OF THE E-HEATING ZONE 

We here predict the location of the e-heating zone in pro¬ 
toplanetary disks using the methodology described in Sec¬ 
tion 3.2. We conduct a parameter study varying the mid¬ 
plane plasma beta Pc, grain size a, dust-to-gas mass ratio fdg, 
and surface density scaling factor f^. Eollowing Sano et al. 
(2000), we select the MMSN (/^ = 1 and q = 3/2) with 
(2 = 0.1 yum, fdg = 0.01, andy^c = 1000 as the fiducial model. 
We start out with this fiducial model in Section 4.1, and dis¬ 
cuss the dependence on the parameters in the subsequent sub¬ 
sections. A summary of the parameter study is given in Ta¬ 
ble 1. We also describe ion heating in Section 4.5. 

4.1. Fiducial Disk Model 

Eigure 3 shows the two-dimensional (radial and vertical) 
map of the dead, active, and e-heating zone in the fiducial 
disk model. The MRI criterion in the ideal MHD limit (Equa¬ 
tion (28)) is satisfied at altitudes below z = Ffideai ~ 23H (see 
Equation (30)). The region above this height is MRI-stable 
with the MRI modes suppressed by too strong magnetic ten¬ 
sion. The dead zone is located inside 24 AU from the cen¬ 
tral star and near the midplane where the gas is shielded from 


ionizing irradiation. The size of the dead zone for this disk 
model is consistent with the prediction by Sano et al. (2000) 
(see their Eigure 7(b)), although their dead zone is slightly 
thicker than ours because of the neglect of X-ray ionization. 
We find that the e-heating zone extends from the outer edge of 
the dead zone out to 82 AU from the central star. This means 
that MRI turbulence can develop without affected by electron 
heating only in the outermost region of r > 80 AU. 

To illustrate how our zoning criteria work in this particular 
example, we plot in Eigure 4 the relation between the current 
density / and electric field E in the midplane at 15 AU, 45 AU 
and 90 AU, which represent the dead, e-heating, and active 
zones, respectively. Recall that for fixed E, MRI turbulence 
grows if J{E) > Ja=i and decays otherwise (Equation (34)). 
At 15 AU, J(E) falls below /a=i for all values of E, imply¬ 
ing that the MRI is unable to grow at this location. At 45 
AU, the MRI growth condition is satisfied during the initial 
growth stage of 10“^^ esu cm“^, but breaks down before 
J reaches /max because of the decrease in J{E) due to elec¬ 
tron heating. This implies that MRI turbulence is allowed to 
grow in the initial stage but saturates at a level lower than that 
for fully developed turbulence. At 90 AU, J{E) reaches /max 
before electron heating sets in, implying that fully developed 
MRI turbulence is sustained here. 

In contrast to electron heating, ion heating is found to be 
negligible at all locations in the fiducial disk model. In the e- 
heating zone, the electric field strength at the saturation point 
is typically < lO^E’cnt (see the center and right panels of Eig¬ 
ure 4), which is an order of magnitude lower than the field 
strength required for ion heating, lO^E’cnt- 

4.2. Dependence on the Magnetic Field Strength 

Eigure 5 shows how the size of the dead and e-heating zones 
depend on the midplane plasma beta Pc. Recall that a higher 
Pc corresponds to a weaker magnetic field B threading the 
disk. As we increase Pc^ the dead zone expands because the 
Elsasser number A oc decreases. On the other hand, we 
find that the boundary between the e-heating and active zones 
is less sensitive to the choice of Pc. As can be inferred from 
the middle and right panels of Eigure 4, this boundary is ap¬ 
proximately determined by the condition that the current den¬ 
sity J{E) reaches /max at a local maximum lying dXE ^ E’crit- 
Since both E’crit and /max are independent of B and hence of 
Pc, so is the boundary between the e-heating and active zones. 

4.3. Dependence on the Grain Size and Dust-to-Gas Mass 

Ratio 

The size and amount of dust grains in disks are important 
parameters in the ionization model as they efficiently remove 
plasma particles from the gas phase. Obviously, these quanti¬ 
ties change as the grains coagulate, settle, or are incorporated 
by even larger solid bodies like planetesimals. We here ex¬ 
plore how the change of these parameters affect the size of 
the dead and e-heating zones. 

To begin with, we show in Eigure 6 the location of the dead, 
active, and e-heating zones with the dust-to-gas ratio fdg fixed 
to 0.01 but with the grain size a varying between 0.1 pm and 
100 pm. We can see that the e-heating zone shrinks with in¬ 
creasing grain size. On increasing a by a. factor of 10, the 
outer radius of the e-heating zone decreases by a factor of 

2. Qualitatively, this is simply because the ionization frac¬ 
tion of the gas increases with decreasing total surface area of 
the grains. Equation (20) shows that the electron abundance 













Figure 4. Relations between the current density J and electric field strength E in the midplane at 15 AU (left panel), 45 AU (middle panel), and 90 AU (right 
panel), which represent the J-E relations in the dead, e-heating, and active zones, respectively. The thick solid line shows the current-field relation /(£), the 
dotted line the maximum current density of MRI, /max (Equation (33)), the vertical gray solid line the criterion for electron heating, Ecrit (Equation (18)), and 
the dashed line the critical current density /a=i below which the MRI decays owing to Ohmic dissipation (Equation (35)). The black dots on the J-E relations 
indicate the saturation points at which either fully developed (/(E) = /max) or self-regulated (/(E) = /a=i) MRI turbulence is sustained. 



Figure 5. Same as Figure 3, but for different values of the midplane beta /5c. 




Figure 6. Same as Figure 3, but for different values of the grain size a. 


Xe = nejun in equilibrium is inversely proportional to the to¬ 
tal surface area of grains per unit volume Ana^rid as long as 
adsorption of plasma particles onto the grains dominate over 
gas-phase recombination. When dust grains aggregate, their 
total surface area decreases inversely proportional to a, and 
hence the electron abundance increases linearly with a. The 
resulting increase in the electric conductivity causes a shift of 
the J-E curve toward higher /, enabling the curve to cross 
the J = /max line at smaller orbital radii. We also find that the 
outer radius of the dead zone decreases at a similar rate to that 


of the e-heating zone when we go from a = fim to 1 yum. 
However, the decrease in the dead zone size stops beyond this 
grain size because gas-phase recombination takes over plasma 
adsorption onto dust grains. As a consequence, the e-heating 
zone becomes narrower and narrower as a increases beyond 
10 yum. 

Decreasing the dust-to-gas mass ratio fdg has a similar ef¬ 
fect to increasing the grain radius because the total surface 
area of the grains is linearly proportional to fdg. This can 
be seen in Figure 7, where we show the location of the dead 
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Figure 7. Same as Figure 3, but for different values of the dust-to-gas mass ratio fdg. 
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Figure 8. Same as Figure 3, but for different values of the surface density scaling factor fz. 


and e-heating zones for <2 = 0.1 yum with fdg varying between 
10“^ and 10“^. We see that the outer radii of the active and e- 
heating zone decrease by a factor of 2 when fdg is decreased 
by a factor of 10. This trend is similar to what we have seen 
when increasing the grain radius by the same factor. 

4.4. Dependence on the Disk Mass 

Finally, we examine how the size of the e-heating zone de¬ 
pends on the disk mass. Figure 8 shows the location of the 
e-heating zone for different values of fz. Here, we fix the 
dust-to-gas mass ratio fz so that both the gas and dust den¬ 
sities scale with fz. We find that the e-heating zone expands 
toward larger orbital radii and higher altitudes as fz increases. 
In the horizontal direction, the expansion is mainly due to the 
increased amount of dust grains with increasing fz. As we 
have explained in 4.3, the ionization fraction of the gas scales 
inversely with Ana^nd, and hence with fz. Therefore, increas¬ 
ing fz by a factor has the same effect as increasing fdg by 
the same factor as long as the ionization rate f is unchanged 
(which is approximately true at ~ 100 AU where cosmic rays 
penetrate down to the midplane). This is exactly what we see 
in Figures 7 and 8, where the e-heating zone expands to 150 
AU when either fdg or fz is increased by the factor of 10 from 
the fiducial value. By contrast, the vertical expansion of the 
e-heating zone is caused by the attenuation of X-rays that oc¬ 
curs at higher altitudes with increasing gas column density. 

4.5. Ion heating 

We observe ion heating in two cases where Pc = 100 and 
where fdg = 0.1. Figure 9 plots the distribution of the ion 




Figure 9. Two-dimensional distribution of the ion temperature Tf for the 
ySc = 100 model (upper panel) and fdg = 0.1 model (lower panel). The solid 
lines show the boundary of the e-heating zone, while the dotted lines show 
'^ideal- 

temperature Tt in the saturated state for these cases. In the 
case of = 100 (the upper panel of Figure 9), U is 3-4 times 
higher than the temperature in a region slightly outside the 
e-heating zone. In this case, the Elsasser number A exceeds 
unity even after electron heating reduces A. This allows the 
electric field strength to reach the critical value for ion heating 
(« 103£eri,) in the vicinity of the e-heating zone. In the case of 
fdg = 0.1 (the lower panel of Figure 9), ion heating takes place 
near the upper boundary of the e-heating zone. However, the 
region is very narrow, and the temperature rise is less than 
2T. Therefore, in this case, ion heating might be practically 
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negligible. 


^max- Assuming for this case, we have 


5. SATURATION OF TURBULENCE IN THE 
E-HEATING ZONE 


6B(i) 


10V2B,0 


J 

Anax 


(39) 


We have shown in Section 4 that self-regulation of the MRI 
due to electron heating can occur over a large region of pro¬ 
toplanetary disks. Then the question arises how strongly the 
e-heating will suppress the MRI turbulence in the e-heating 
zones. This question can only be fully addressed with MHD 
simulations including magnetic diffusion and electron heating 
in a self-consistent manner, which is far beyond the scope of 
this study. In this section, we attempt to estimate the satura¬ 
tion level of MRI turbulence from simple scaling arguments. 

As usual, we quantify the strength of turbulence with the 
Shakura-Sunyaev a parameter a = TrcfylP, where P = pc^ 
is the gas pressure and is the component of turbulent 
stress. In MRI-driven turbulence, is generally dominated 
by the turbulent Maxwell stress -SBrSBfp/An (Hawley et al. 
1995; Miller & Stone 2000), where SBr and SB^p are the radial 
and azimuthal components of the turbulent (fluctuating) mag¬ 
netic fields. Therefore, we evaluate the a parameter for MRI 
turbulence as 


O^MRl 


SB^SBfp 

4npc^ 


(36) 


In reality, the Reynolds stress (Fleming & Stone 2003; 
Okuzumi & Hirose 2011) or the coherent component of the 
Maxwell stress (e.g.. Turner & Sano 2008; Gressel et al. 2011) 
can dominate over the turbulent Maxwell stress at locations 
where the MRI is significantly suppressed. However, we do 
not include these components in our c^mri because they do not 
reflect the local MRI activity at such locations (see the refer¬ 
ences above). 

Next we relate the amplitude of turbulent magnetic fields to 
the amplitude of the electric current density J = \J\ using the 
Ampere’s law J = (c/4;r)V x SB. We neglect large-scale, co¬ 
herent components in B since the electric current is inversely 
proportional to the length scale of fields. We assume that 
the magnetic field in MRI-driven turbulence is dominated by 
the azimuthal component SB(p and varies over a length scale 
~ dideah where dideai is the wavelength of the most unstable 
MRI modes already introduced in Equation (29). Then, from 
the Ampere’s law, one can estimate the magnitude of the cur¬ 
rent density as 


/=—IVx^l 

An 


An Q. ^ 



(37) 


where we have replaced the derivative V with wavenumber 
27r//lideai = ^/^Az- If we use the maximum current /max 
for fully developed MRI turbulence (Equation (33)), Equa¬ 
tion (37) results in a simple scaling relation 


SB^ r- J 

. (38) 


For fully developed MRI turbulence where J ^ /max, the 
above equation predicts SB(plBz ~ 10, in agreement with the 
results of MHD simulations (e.g., Hawley et al. 1995; Sano 
et al. 2004). 

Now let us consider situations where e-heating is so effec¬ 
tive that the growth of the MRI is saturated at / ^^ /a=i 


This equation predicts the amplitude of SB(p as a function 
of Bzo and ///max- MHD simulations show that SBr ~ 
-(O.A.. .0.6)SB(p in MRI turbulence (Hawley et al. 1995; 
Sano et al. 2004). Assuming that this scaling also holds in 
our case, we have SBySBfp ^ -1005^Q(///max)^. Finally, sub¬ 
stituting this into Equation (36), we obtain the scaling relation 
between aum and ///max, 


auRi ~ 


100^2 

Anpc^ 



^0.2 




(40) 


where JSq = ^npc^JB^^ = Pc exp {-z^l2H^) is the plasma beta 
(not necessarily at the midplane) associated with the net verti¬ 
cal field Bzo. Formally, the derivation leading to Equation (40) 
breaks down when MRI is so active that SB^ » B^q and 
/ ~ /max- Nevertheless, we find that Equation (40) reproduces 
the results of ideal MHD simulations with a reasonably good 
accuracy. Equation (40) predicts that ckmri ~ 2 for JSq = 10^ 
and auRi ~ 0.02 for JSq = 10^ when / = /max- These are 
consistent with the results of isothermal simulations by Sano 
et al. (2004) showing that the Maxwell component of a is ~ 1 
for JSo = 10^ and ~ 0.01 for JSq = 10"^ (see their Table 2, col¬ 
umn (10)). Therefore, we will apply Equation (40) to both the 
e-heating zone and active zone. 

The left panel of Figure 10 show the radial distribution of 
o^MRi for the fiducial disk model predicted from Equation (40). 
Here we plot the midplane value (TMRi,mid = o^mrKz = 0) and 
the density-weighted average in the vertical direction, 

Xh'""" o;MRi{z')p{z')dz' 

auKi = ^-, ( 41 ) 

where we have assumed auRi = 0 in the magnetically dom¬ 
inated atmosphere at |z| > /^ideai- The former quantity mea¬ 
sures the MRI activity at the disk midplane, while the latter 
quantity is more closely related to the vertically integrated 
mass accretion rate (Suzuki et al. 2010). For the fiducial 
disk model, we find that aMRi,mid ~ 10 ^ and 10 ^ at the 
inner and outer edge of the e-heating zone (20 AU and 80 
AU), respectively. These values are more than two orders 
of magnitude lower than the value aMRi,mid = 0.2 in the ac¬ 
tive zone (r > 80 AU). This implies that the MRI is “virtu¬ 
ally dead” deep inside the e-heating zone. We also find that 
^rMRi,mid changes discontinuously at the boundary between the 
e-heating and active zones. The reason is that when the sat¬ 
urated state changes at the point, ///max also changes from 
unity to one order of magnitude because of the N-shaped 
current-field relation (see middle and right panels of Figure 
4). The vertical average q^mri decreases more slowly with de¬ 
creasing r, because the upper layer of the disk remains MRI- 
active (see Figure 3). This picture is qualitatively similar to 
the classical layered accretion model of Gammie (1996). In 
right panel of Figure 10, we also plot the radial distribution 
of aMRi,mid and auRi for a disk with Pc = 10"^. We find that 
o^MRXmid in e-heating zone is almost unchanged from the fidu¬ 
cial disk. The reason is that increase of (///max)^ ~ 10 cancels 
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Figure 10. Radial distribution of c^mri (Equation (40)) for the fiducial model (left panel) and fie = 10^ (right panel). The solid black line shows aMRi,mid 
including electron heating on the mid-plane, and the solid blue line shows umri including electron heating integrated in the z-direction. The dashed black line 
shows crMRi,mid without including electron heating on the mid-plane, and the dashed blue line shows q^mri without including electron heating integrated in the 
z-direction. 


out the depletion of ^ 10“Mn Equation (40). Therefore, 
<TMRi,mid remains low saturation level. 

In summary, our simple estimate predicts that MRI turbu¬ 
lence can be significantly suppressed in the e-heating zone. In 
this sense, the e-heating zone acts as a extended dead zone. 
However, our estimate relies on the hypothetical scaling be¬ 
tween the and turbulent Maxwell stress and ///max, which is 
as yet justified by MHD simulations.^ In order to test our pre¬ 
diction, we will perform resistive MHD simulations including 
electron heating in future work. 

6. CHARGE BARRIER AGAINST DUST GROWTH IN 
THE E-HEATING ZONE 

So far we have focused on the role of electron heating on the 
saturation of MRI turbulence. As pointed out by 0115, elec¬ 
tron heating also has an important effect on the growth of dust 
grains. In an ionized gas, dust grains tend to be negatively 
charged because electrons collide and stick to dust grains 
more frequently than ions. The resulting Coulomb repul¬ 
sion slows down the coagulation of the grains through Brow¬ 
nian (thermal) motion. This “charge barrier” is also present 
in weakly ionized protoplanetary disks, in which dust grains 
tend to be charged as in a fully ionized gas when their size 
is larger than 1 fim (Okuzumi 2009; Matthews et al. 2012). 
The important role of electron heating in this context is that 
heating electrons further promote the negative charging of the 
grains, because the grain charge in a plasma is linearly pro¬ 
portional to the electron temperature (e.g., Shukla & Mamun 
2002). In this section, we explore how this affects dust coag¬ 
ulation in the e-heating zone. 

Eor simplicity, let us assume that dust grains have the single 
radius a and charge Z. The grains can collide with each other 

^ However, there are some support for Equation (40) from MHD simu¬ 
lations including ambipolar diffusion, not Ohmic dissipation. Bai & Stone 
(2011) reported the Maxwell component of a (their Table 2) and the cu¬ 
mulative probability distribution of J (Figure 6) for three simulation runs 
with pQ = 400 and with different values of ambipolar diffusivity. Their 
results show that UMaxweii ~ 0-17, 0.029, and 0.0041 for models with 
i/imay ~ 1,0.3, and 0.1 (median values), respectively. These are consistent 
with Equation (40) predicting that ctmri ~ 0.5, 0.045, and 0.005 for these 
values of // /max • 


if the condition 


^col ^ ^elc 


(42) 


is satisfied (Okuzumi 2009). Here, 6coi is the kinetic energy 
of the relative motion of two colliding grains, and 


^elc 


{eZf 

2a 


(43) 


is the Coulomb repulsion energy of the grains just before con¬ 
tact. We focus on small dust grains near the midplane and as¬ 
sume that the relative motion is dominated by Brownian mo¬ 
tion and turbulence-induced motion. Then, the kinetic energy 
of relative motion can be expressed as 


^col ~ ^Brown ^turb, 


(44) 


where ^Brown and 6turb are the kinetic energy of Brownian mo¬ 
tion and turbulence-induced motion, respectively. Brownian 
motion is the thermal motion of grains, and ^Brown is approx¬ 
imately expressed as 


^Brown ~ 2^^th’ 


(45) 


where the thermal velocity of grains Wth is expressed as Wth = 
VSkr !Tim and the reduced mass of grains ja is expressed as 
II = m^/(m m) = mil. The relative energy of turbulence- 
induced motion is expressed as 


i 2 

^turb ~ 2^^^^turb) , (46) 

where Ai/turb is the relative velocity of the grains excited by 
turbulence. Eor small grains, Ai/turb is approximately given by 
(Weidenschilling 1984; Ormel & Cuzzi 2007) 

^f^turb ~ V^dispR^ ^ Cs^Ts, (47) 


where Ofdisp = /cj is the velocity dispersion of the gas 

normalized by cj. Re is the Reynolds number of turbu¬ 
lence, and 


^ = p.ali ^j^Csp) 


T. 


(48) 
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Figure 11. Effectiveness of the charge barrier against grain growth as a 
function of the grain size at the midplane 35 AU in the fiducial model. The 
solid line (red) shows Seic/^coi including electron heating, and the dashed 
line (blue) shows 6eic/^coi without including electron heating. The horizon¬ 
tal dotted line shows 6eic/^coii = above which a strong Coulomb repul¬ 
sion between the grains suppresses their mutual collision cross section. Here 
it is assumed that adisp(= equal to ctmrl the normalized local 

Maxwell stress given by Equation (36) (but see also Eigure 12). 

is the stopping time of the grains (we have adopted Epstein’s 
drag law for r^). The Reynolds number is expressed as 
Re = ffdispt^?f^~VTmoh where Vmoi is the molecular viscos¬ 
ity. We estimate a with and without electron heating, using 
Equation (40) presented in Section 5. Turbulence dominates 
the collisional energy when is high and/or a is large. Eor 
the moment, we simply assume cfdisp = where auRi 

is the normalized local Maxwell stress introduced in Equa¬ 
tion (36). This assumption holds when the Reynolds stress in 
the e-heating zone is comparable to the Maxwell stress. In re¬ 
ality, the Reynolds stress in the e-heating zone might be higher 
than the Maxwell stress for a reason discussed later. There¬ 
fore, the estimate of 6turb presented here should be taken as a 
lower limit. 

To obtain Z and c^mri, we calculate the ionization fraction 
(Section 2.3), determine the turbulent state (Section 3.2), and 
estimate the MRI-turbulent viscosity (Section 5) with chang¬ 
ing grain radius a at a location. We then obtain 6coi and 6eic 
by above-mentioned method. It should be noted that grains 
have single size and changing grain radius means changing 
the size of all grains at the location. Thus, the turbulent state 
at the location also depends on a. 

In Eigure 11, we plot the ratio 6eic/^coi as a function of 
(2 at 35 AU in the midplane for the fiducial disk model. The 
ratio quantifies the effectiveness of the charge barrier: the col¬ 
lisional cross section of two equally charged grains is sig¬ 
nificantly suppressed when Sq\c/Sco\ » 1. We find that 
electron heating significantly enhances the charge barrier for 
submicron-sized grains. If electron heating is not included, 
this location belong to the dead zone and the active zone with 
grain size being <0.05 /im and >0.05 /im, respectively. In 
this case, 6eic/^coi is much lower than unity in all a. Thus we 
can conclude that dust grains at this location can grow with¬ 
out the charge barrier. On the other hand, if electron heating 
is included, this location belongs to the e-heating zone when 


Figure 12. Same as Figure 11, but we here evaluate Udisp = using 

Equation (49) assuming that sound waves propagate from upper MRI active 
layers to the midplane as it is observed for the conventional dead zone. 

0.05 yum < a < 1.4 fim (see also Eigure 6). In the e-heating 
zone, grains are charged by heated electrons, leading to in¬ 
crease of 5eic, and MRI turbulence as collisional source is 
well suppressed, leading to decrease of 6turb- Consequently, 
^eic/^coi is larger than unity when 0.08 /im < a < 0.5 /im. 
In particular, 5eic/^coi takes its maximum value of 40 at 
a = 0.2 /im corresponding to ^Brown = ^turb- Both the sup¬ 
pression of turbulence and grain charge would enhance the 
charge barrier. 

There are at least two mechanisms that could drive further 
growth of dust in the e-heating zone. One is vertical turbulent 
mixing of dust particles as already pointed out by Okuzumi 
et al. (2011). In general, the charge barrier is less significant 
at higher altitudes where dust particles have a higher colli¬ 
sion energy due to vertical settling (and due to if MRI is ac¬ 
tive there). Electron heating, which was not considered by 
Okuzumi et al. (2011), does not change this picture because it 
is also ineffective at high altitudes. Micron-sized grains in the 
e-heating zone can easily be lifted up to such high altitudes if 
only weak turbulence is present there (Turner et al. 2010, see 
also dust scale height Hd in Section 7.1). The lifted grains are 
allowed to collide and grow there until they fall back to the 
e-heating zone. In this way, small grains in the e-heating zone 
are able to continue growing on a timescale much longer than 
vertical diffusion timescale. Okuzumi et al. (2011) showed 
that the charge barrier is overcome on a timescale of 10^-10^ 
yr, but they did not consider the amplification of grain charg¬ 
ing due to electron heating. How much the growth is delayed 
in the presence of electron heating should be studied in future 
work. 

Another potentially important mechanism is dust stirring by 
random sound waves. It is known that the Reynolds stress 
in a dead zone exceeds the Maxwell stress because of sound 
waves propagating from upper MRI-active layers (e.g., Ero- 
mang & Papaloizou 2006; Turner et al. 2010; Okuzumi & Hi- 
rose 2011). If this is also the case for our e-heating zone, the 
assumption a^disp = <tmri would significantly underestimate 
the particle collision energy in the e-heating zone. To estimate 
this effect, we now calculate a^disp using an empirical formula 
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for the gas velocity dispersion in the dead zone (Okuzumi & 
Hirose 2011), 


{Sv^) ^ 0.78aMRic? exp j, (49) 

where q^mri is the density-weighted vertical average of q^mri 
defined by Equation (41). Equation (49) expresses the am¬ 
plitude of random sound waves inside a dead zone. Eigure 
12 shows 6eic/^coi in this case and is obtained in the same 
way as in Eigure 11 but we here use Equation (49) for c^disp in 
Equation (47). The use of Equation (47) for the sound wave- 
driven collision velocity assumes that the time correlation of 
the waves’ velocity fluctuations exponentially decays on the 
timescale of as in the Kolmogorov turbulence. We find 
that SqIcIScoI now falls below unity at all grain sizes. Thus, 
sound waves traveling from MRI-active layers, if they exist, 
could help dust overcome the charge barrier in the e-heating 
zone. However, the argument made here is not conclusive be¬ 
cause the induced collision velocity depends on the assumed 
time correlation function, or equivalently power spectrum, of 
the random sound waves. If the power spectrum of the waves 
has only a small amplitude at high frequencies (to which small 
dust particles are sensitive) compared to the turbulent spec¬ 
trum, the wave-induced collision velocity would be lower than 
given by Equation (47). The spectrum of velocity fiuctuations 
in the e-heating zone should be studied in future MHD simu¬ 
lations. 


7. DISCUSSION 
7.1. Dust Diffusion 

We have assumed so far that the dust-to-gas mass ratio is 
vertically constant. This assumption breaks down when dust 
particles settle toward the midplane. If this is the case, the 
dust-to-gas ratio would decrease at high altitudes, and conse¬ 
quently the e-heating zone would shrink in the vertical direc¬ 
tion as expected from Eigure 7. 

However, as we will show below, dust settling is negligible 
even in the e-heating zone because even weak turbulence is 
able to diffuse small grains to high altitudes. Youdin & Lith- 
wick (2007) analytically derived dust scale height in the 
sedimentation-diffusion equilibrium. If the particle stopping 
time Ts is much smaller than the Keplerian timescale 
which is true for small particles, the dust scale height can be 
approximately written as 


/ St 

- , (50) 

\ <^disp,z / 

where St = is the so-called Stokes number and adisp,z = 
I is the vertical component of the velocity dispersion 

normalized by cj. Equation (50) implies that dust settling 
takes place {Hd < H) when St > adisp,z- Under the disk model 
employed in this study, St can be expressed as 


St = 3 X 10"^ 



Therefore, for <2 = 0.1 yum, dust settling in the e-heating zone 
(r ~ 10-100 AU) occurs only if afdisp,^ < 10“^-10“^. In the 
e-heating zone, c^mri ~ 10“^-10“^ at the midplane (see Eig¬ 
ure 10), and therefore we may safely neglect dust settling 
even if the Reynolds stress is as small as the Maxwell stress 


(Q^disp,z ~ o^MRi). A larger a does not change this conclusion, 
because we then would have a higher ckmri or the e-heating 
zone would vanish. 


7.2. Effects of Grain Size Distribution and Porosity 

We have characterized dust grains with a single particle size 
a assuming that the size distribution of dust grains is narrow. 
Under this assumption, the e-heating zone covers only a small 
part of protoplanetary disks when the particles grow to mil¬ 
limeter sizes (see Eigure 6). However, caution is required 
in applying our results to more general cases where particles 
have a size distribution. In such cases, the smallest grains tend 
to dominate the total surface area of dust (which controls the 
ionization balance), whereas the largest grains tend to domi¬ 
nate the total mass of dust, simply because smaller grains have 
a larger area-to-mass ratio. Therefore, it is not obvious what 
the typical particle size is in these cases. 

Here we discuss more quantitatively how we can apply the 
results of single-size calculations to cases with a size distribu¬ 
tion. Let us assume that the particle size distribution is given 
by the power-law form 


dnd _ ^pfdg 
da 


(52) 


with <2min < a < <2max (<^min ^max), where dnd/da is the 
number density of dust particles per unit particle radius, and 
<2min and <2max are the minimum and maximum particle sizes, 
respectively. The distribution is normalized so that that the 
total particle mass density f md(dnd/da)da becomes equal to 
pfdg. Equation (52) applies when the particle size distribution 
is determined by fragmentation cascade (Dohnanyi 1969) and 
is also known to reproduce the size distribution of interstellar 
dust grains (Mathis et al. 1977). The quantity we are inter¬ 
ested in is the total surface area of the particles as it mainly 
determines the ionization balance in a gas-dust mixture (e.g., 
Sano et al. 2000). This can be calculated as 



4na^^^da 

da 


^Pfdg 1 
P* s/amin^max 


(53) 


Note that the factor 1 / sja^m comes from the fact that the in¬ 
tegration in Equation (53) is dominated by the smallest par¬ 
ticles (because a^(dnd/da)da oc d(a~^-^)), whereas the factor 
1 / V^max from the fact that the total mass is dominated by 
the largest particles. By contrast, if all dust particles have a 
single size <2singie, their total surface area is ^na^^^^^^nd,single = 
^pfdg/ip•(^single)- Comparing this with Equation (53), we find 
that the total surface area of particles whose size distribution 
is given by Equation (52) is equal to that of single-size parti¬ 
cles if 


^single ~ 


min‘^max- 


(54) 


Since the total surface area approximately determines the ion¬ 
ization state. Equation (54) may be used to generalize the re¬ 
sults presented in this study to the cases where the particle 
size distribution obeys Equation (52). 

Observations of millimeter dust emission from protoplan¬ 
etary disks suggest that the largest dust particles in the disk 
have a size of centimeters (e.g., Testi et al. 2003; Natta et al. 
2004; Rodmann et al. 2006; Ricci et al. 2010). Assuming 
<^max = 1 cm and <2niin = 0.1 yum, we obtain <2singie = 30 pm. 
In this case, we expect from Table 1 that the e-heating zone 
extends to ~ 15 AU. Thus, even if cm-sized grains exist in 
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protoplanetary disks and the total mass of grains is dominated 
by such large grains, the e-heating zone can be present in the 
disks. 

For the same reason, large dust particles can alone provide 
a large e-heating zone if the dust particles are highly fluffy 
aggregates of tiny grains. Okuzumi (2009) showed that the 
ionization balance is insensitive to the particle radius when the 
fractal dimension is 2, for which the total surface area of the 
aggregates is approximately conserved during the aggregation 
process. 

7.3. Hall Effect and Ambipolar Diffusion 

The plasma heating model employed in this study neglects 
the effects of magnetic fields on the motion of plasma par¬ 
ticles. In terms of non-ideal magnetohydrodynamics, this is 
equivalent to neglecting ambipolar diffusion and Hall effect 
(see, e.g., Wardle 1999). A full treatment of these non-Ohmic 
effects introduces to the model additional complexities aris¬ 
ing from the relative angle between the magnetic and electric 
fields (Okuzumi, Mori, & Inutuska, in prep.), which is beyond 
the scope of this paper. In this subsection, we only briefly dis¬ 
cuss how plasma heating and these non-ideal MHD effects 
could affect each other. 

Ambipolar diffusion can suppress MRI in low density re¬ 
gions of protoplanetary disks (e.g., Blaes & Balbus 1994; 
Hawley & Stone 1998; Kunz & Balbus 2004; Desch 2004; 
Bai & Stone 2011; Simon et al. 2013a,b). If MRI is effectively 
suppressed in the e-heating zone, electric fields may not suf¬ 
ficiently grow to cause electron heating. The effectiveness of 
ambipolar diffusion is characterized by the ambipolar Elsasser 
number Am = ytPilD. (e.g., Blaes & Balbus 1994; Lesur et al. 
2014), where yi = {cTinVin) /(nin + Mi) and pt = mini. Ac¬ 
cording to MHD simulations including ambipolar diffusion, 
MRI-driven turbulence behaves as in the ideal MHD limit if 
Am » 1, while ambipolar diffusion suppresses turbulence if 
Am 1 (e.g., Bai & Stone 2011). Table 2 lists the values 
of Am as well as the ion abundance xi = nt/nn at the inner 
and outer edges of the e-heating zone before electron heating 
sets in (E = 0). We find that Am 0.2-0.7, implying that 
ambipolar diffusion would moderately affect MRI turbulence 
in the e-heating zone. Therefore, MHD simulations includ¬ 
ing both electron heating and ambipolar diffusion are needed 
to assess which effect determines the saturation amplitude of 
MRI turbulence in these outer regions of the disks. 

The Hall effect is also important at r ~ 10-50 AU (see 
Figure 1 of Turner et al. 2014). The Hall effect can either 
damp or amplify magnetic fields, which depends on the rela¬ 
tive orientation between the disk’s magnetic field and rotation 
axis and on the sign of the Hall conductivity (e.g., Bai 2014; 
Wardle & Salmeron 2012). At relatively high gas densities 
(jin > 10^^ cm“^), the Hall conductivity is usually positive 
(Wardle & Ng 1999; Nakano et al. 2002; Salmeron & Wardle 
2003), but can become negative when the number density of 
electrons is significantly lower than that of ions. Interestingly, 
our preliminary investigation shows that the Hall conductivity 
can indeed become negative as the electron number density is 
decreased by electron heating (Okuzumi et al., in prep.). This 
suggests that electron heating might reverse the role of the 
Hall term. Whether this occurs under conditions relevant to 
protoplanetary disks will be studied in future work. 

8. SUMMARY 

We have investigated where in protoplanetary disks the 
electron heating by MRI-induced electric fields affects MRI 


turbulence. Our previous study (0115) showed that elec¬ 
tron heating causes a reduction of the electron abundance, 
and hence an amplification of Ohmic dissipation, when the 
recombination of plasma mainly takes place on dust grains 
rather than in the gas phase. To study where in disks this 
effect becomes important, we constructed a simplified ion¬ 
ization model that takes into account both recombination on 
dust grains and electron heating. The presented model is 
computationally much less expensive than the original elec¬ 
tron heating model by 0115 and allows us to study the ef¬ 
fects of electron heating for a wide range of model parame¬ 
ters. We then searched for locations in a disk where the en¬ 
hanced Ohmic diffusivity limits the saturation level of MRI 
turbulence, which we call the “e-heating zone,” by using ana¬ 
lytic criteria for MRI growth. Our results can be summarized 
as follows: 

1. We find that the e-heating zone can cover a large part of 
a protoplanetary disk when tiny dust grains are abun¬ 
dant. For instance, in a minimum-mass solar nebula 
with 1% of its mass consisting of 0.1-yum-sized dust 
grains, the e-heating zone extends out to 80 AU from 
the central star (Figure 3; Section 4.1). In this case, 
MRI turbulence can develop without being affected by 
electron heating only in the outermost region of r > 
80 AU. 

2. In the e-heating zone, the saturation level of MRI tur¬ 
bulence is expected to be considerably lower than that 
in fully MRI-active zones because the electron heating 
sets an upper limit to the electric current density attain¬ 
able in MRI turbulence. Our simple estimate based on 
scaling arguments (Section 5) predicts that for our fidu¬ 
cial disk model, the turbulence a parameter for MRI 
turbulence should be reduced to ~ 10“^ and 10“^ at 
the inner and outer edges of the e-heating zone, respec¬ 
tively (Figure 10). This implies that the MRI is “virtu¬ 
ally dead” deep inside the e-heating zone. 

3. Dust grains in the e-heating zone acquire a high nega¬ 
tive charge due to the frequent collisions with electri¬ 
cally heated electrons. This strengthen the charge bar¬ 
rier against the growth of micron-sized grains originally 
predicted by Okuzumi (2009) (Figure 11; Section 6). 
At midplane 35 AU in the fiducial model, the electric 
repulsion energy is larger than the collisional energy 
when the grain size is in the range of ~ 0.08-0.5 pm. 
We find that electron heating significantly enhances the 
charge barrier for submicron-sized grains. 

Our estimate of the turbulence strength in the e-heating 
zone (Equation (40)) largely relies on the scaling relations 
between turbulent quantities observed in previous MHD sim¬ 
ulations. Although these scalings well predict the saturation 
level of MRI turbulence without electron heating, it is unclear 
whether they are still valid even in the presence of electron 
heating. Our future work will address this issue by perform¬ 
ing MHD simulations including electron heating. We have 
also neglected the effect of magnetic fields on the kinetics of 
plasma, which means that non-Ohmic effects such as the Hall 
effect and ambipolar diffusion are excluded from our analysis. 
However, these effects generally overwhelm Ohmic diffusion 
in outer parts of protoplanetary disks. Our estimate indicates 
that ambipolar diffusion would moderately suppress MRI in 
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Table 2 

Am and Ion Abundance x/ in E-heating Zone for Various Parameter Sets 



a ( jdm ) 

fdg 

h 

Am in e-heating zone 
Inner edge Outer edge 

Xi in e-heating zone 
Inner edge Outer edge 

10 ^ 

0.1 

10-2 

1 

0.14 

0.56 

1 . 9 x 10-12 

4 . 7 x 10 -" 

10 ^ 

0.1 

10-2 

1 

0.17 

0.62 

3.4x10-'2 

5 . 9 x 10 -" 

10 “* 

0.1 

10-2 

1 

0.23 

0.62 

7.4x10-‘2 

5 . 9 x 10 -" 

105 

0.1 

10-2 

1 

0.41 

0.62 

2 . 5 x 10 -" 

5 . 9 x 10 -" 

10 ^ 

0.1 

10-2 

1 

0.17 

0.62 

3 . 4 x 10 -'^ 

5 . 9 x 10 -" 

10 ^ 

1 

10-2 

1 

0.21 

0.72 

1.7x10-‘2 

2 . 8 x 10 -" 

10 ^ 

10 

10-2 

1 

0.43 

0.82 

2.3x10-'2 

1 . 3 x 10 -" 

10 ^ 

100 

10-2 

1 

0.54 

0.72 

2.6x10-'2 

5.6x10-‘2 

10 ^ 

0.1 

10-1 

1 

0.16 

0.57 

8.5x10-'2 

1 . 2 x 10 -'“ 

10 ^ 

0.1 

10-2 

1 

0.17 

0.62 

3.4x10-‘2 

5 . 9 x 10 -" 

10 ^ 

0.1 

10-3 

1 

0.24 

0.71 

2.1x10-‘2 

2 . 8 x 10 -" 

10 ^ 

0.1 

10-4 

1 

0.41 

0.80 

2.3x10-'2 

1 . 3 x 10 -" 

10 ^ 

0.1 

10-2 

10 

0.32 

1.26 

1.8x10-'2 

2 . 5 x 10 -" 

10 ^ 

0.1 

10-2 

3 

0.22 

0.82 

2.5x10-‘2 

3 . 9 x 10 -" 

10 ^ 

0.1 

10-2 

1 

0.17 

0.62 

3.4x10-'2 

5 . 9 x 10 -" 

10 ^ 

0.1 

10-2 

0.3 

0.16 

0.65 

5.4x10-'2 

9 . 5 x 10 -" 


the e-heating zone (Am 0.2-0.7; Section 7.3). Therefore, 
MHD simulations including both electron heating and non- 
Ohmic resistivities will be needed to assess which effect de¬ 
termines the saturation amplitude of MRI turbulence in outer 
regions of the disks. We will address these open questions 
step by step in future work. 
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